Analogs of quantum Hall effect edge states in photonic crystals 
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"Photonic crystals" built with time-reversal-symmetry-breaking Faraday-effect media can exhibit 
"chiral" edge modes that propagate unidirectionally along boundaries across which the Faraday 
axis reverses. These modes are precise analogs of the electronic edge states of quantum Hall effect 
(QHE) systems, and are also immune to backscattering and localization by disorder. The "Berry 
curvature" of the photonic bands plays a role analogous to that of the magnetic field in the QHE. 
Explicit calculations demonstrating the existence of such unidirectionally-propagating photonic edge 
states are presented. 



I. INTRODUCTION 

The control of the flow of light using photonic band-gap 
(PBG) materials has received considerable attention over 
the past decade 1 . Moreover, the potential for using arti- 
ficially structured "metamaterials" , such as the recently 
discovered "left-handed media"—, has shown considerable 
technological promise. In the past, significant progress 
has been achieved in the field of photonics by making use 
of analogies with electronic systems. For instance, the 
typical PBG material, a system with a spatially vary- 
ing and periodic dielectric function, was motivated by 
the well known physics of electronic Bloch states; the 
dielectric scattering of light in periodic media presents 
the same formal solutions as those for the scattering of 
electrons in periodic potentials. 

Previous photonic bandstructure calculations have fo- 
cused on the frequency dispersion of the photon bands; 
it has been usually been assumed that a knowledge of 
the spectrum alone represents a complete understanding 
of dynamics of the system. A primary goal of such cal- 
culations has been the quest for a PBG material with a 
complete bandgap throughout the Brillouin zone in some 
frequency range, which would prevent the transmission 
of light with frequency in the range of the Band gap. 
Both two and three dimensional bandstructures possess- 
ing these properties have now been discovered^. 

Recently, however, in the study of electronic systems, it 
has become apparent that, even in the absence of interac- 
tion effects, the dispersion relations of the energy bands 
do not fully characterize the semiclassical dynamics of 
wavepackets, unless both spatial inversion symmetry and 
time-reversal symmetry are unbroken . The additional 
information, which is not obtainable from knowledge of 
the energy bands e n (fc) alone, is the variation of the Berry 
curvature^ J-^ b (k) — e ahc f2 nc (fc), which is an antisym- 
metric tensor in fc-space, where O rl (fc) is analogous to a 
"magnetic field" (flux density) in fc-space. The "Berry 
curvature" in fc-space is related to the Berry phased in 
the same way that the Bohm-Aharonov phase of an elec- 
tronic wavepacket is related to the magnetic flux density 
in real space. 

While the uniform propagation of wavepackets in per- 
fectly translationally-invariant systems does not involve 



the Berry curvature, the "semiclassical" description of 
the acceleration of wavepackets in media with spatial in- 
homogeneity of lengthscales large compared to the under- 
lying lattice spacing is incomplete if it is not taken into 
account. Recently Onoda et al£ have pointed out the role 
of Berry curvature in photonic crystals without inver- 
sion symmetry; while these authors characterize this as a 
"Hall effect of light" , the Hall effect in electronic systems 
is associated with broken time-reversal symmetry rather 
than broken spatial inversion symmetry, and we have re- 
cently discussed^ some of the at-first-sight-surprising ef- 
fects that broken time-reversal symmetry could produce 
in photonic systems. 

In the presence of non-vanishing Berry curvature, the 
usual "semiclassical" expression for the group velocity 
of the wavepacket is supplemented by an additional 
"anomalous" contribution proportional to its accelera- 
tion and the local Berry curvature of the Bloch band. 
(The semiclassical treatment of electron dynamics be- 
comes ray optics in the photonic context). This "anoma- 
lous velocity" has played an important role in under- 
standing recent experiments on the anomalous Hall effect 
of ferromagnets^, for example. 

Perhaps the most remarkable among the "exotic" ef- 
fects associated with Berry curvature, however, is the 
quantum Hall effect, which has been the focus of in- 
tensive experimental and theoretical study in condensed 
matter physics for over two decades. The physics of the 
quantum Hall regime and its connection with Berry Cur- 
vature phenomena is now well understood. The possibil- 
ity of transcribing some of the main features of the quan- 
tum Hall effect to photonic systems, which brings into 
play new possibilities in photonics, is the topic of this pa- 
per. Specifically, we shall concern ourselves with analogs 
of "chiral" (unidirectional) quantum Hall edge states in 
photonic systems with broken time-reversal symmetry. 

The quantum Hall effect is usually associated with two 
dimensional electron systems in semiconductor hetero- 
junctions in strong applied magnetic fields. By treat- 
ing the plane of the heterojunction as a featureless two- 
dimensional (2D) continuum, and considering the quan- 
tum mechanical motion of electrons in the presence of a 
magnetic field, one obtains the electronic Landau levels. 
The key feature giving rise to the quantization of the Hall 
conductance is the incompressibility of the electron fluid: 



2 



either due to the Pauli principle at integer Landau level 
fillings, with the spectral gap to the next empty level 
given by the cyclotron frequency, or when a gap opens 
due to strong electron-electron interactions at fractional 
fillingsii. 

While in the experimentally-realized systems, the 
quantum Hall effect derives from a strong uniform com- 
ponent of magnetic flux density normal to the 2D plane 
in which the electrons move, the integer quantum Hall 
effect can also in principle derive from the interplay of a 
periodic bandstructure with a magnetic field. Externally 
applied in periodic structures give rise to the celebrated 
Hofstadter model of Bloch bands with an elegant fractal 
spectral structure depending on the rational value of the 
magnetic flux through the unit celii2A£. The influence of 
the lattice on the quantum Hall effect was further investi- 
gated in an important paper by Thouless et al (TKNN) 14 , 
who discovered a topological invariant of 2D bandstruc- 
tures known as the "Chern Number" a quantity that was 
later interpreted in terms of Berry curvature^. 

At first sight, it seems implausible that any of the phe- 
nomena associated with the quantum Hall effect can be 
transcribed to photonics, fncompressibility and Landau 
level quantization require fermions and charged particles, 
respectively, and it is not clear how one could introduce 
an effect similar to the Lorentz force due to a magnetic 
field on a system of neutral bosons. However, a hint that 
possible analogs could exist in photonics, comes from the 
fact that a "zero-field" quantum Hall effect without any 
net magnetic flux density (and hence without Landau lev- 
els) could occur in systems consisting of "simple" Bloch 
states with Broken time-reversal symmetry, as was shown 
some time ago by one of U9~. The explicit "graphene- 
like" model investigated in Ref. ^| exploits the topo- 
logical properties of Bloch states, which motivated us to 
construct its photonic counterpart. This model has also 
turned out to be a very useful for modeling the anoma- 
lous Hall effect in ferromagnetic metals^, and a recently 
proposed "quantum spin Hall effect"—. 

While incompressibility of the fluid in the bulk quan- 
tizes the Hall conductance, perhaps the most remarkable 
feature of quantum Hall systems is the presence of zero 
energy excitations along the edge of a finite system. In 
these edge states, electrons travel along a single direction: 
this "one-way" propagation is a symptom of broken time- 
reversal symmetry. In the case that the integer quantum 
Hall effect is exhibited by Bloch electrons, as in the Hof- 
stadter problem studied by TKNN 14 , it is related to the 
topological Chern invariant of the one-particle bands and 
therefore. The edge states necessarily occur at the inter- 
face between bulk regions in which there is a gap at the 
Fermi energy, which have different values of the sum of 
the Chern invariants of the fully occupied bands below 
the Fermi level. While the integer quantum Hall effect 
in such a system itself involves the filling of these bands 
according to the Pauli principle, and hence is essentially 
fermionic in nature, the existence of the edge states is 
a property of the one-electron band structure, without 



reference to the Pauli principle, which suggests that this 
feature is not restricted to fermionic systems. We have 
found that they indeed have a direct photonic counter- 
part, and this leads to the idea that topologically non- 
trivial unidirectionally propagating photon modes can oc- 
cur along domain walls separating two PBG regions hav- 
ing different Chern invariants of bands below the band 
gap frequency. In this paper, we present the formal ba- 
sis of such modes, along with explicit numerical exam- 
ples, simple model Hamiltonians, and semiclassical cal- 
culations confirming the concept. 

We note finally, that while Berry phases are usually as- 
sociated with quantum mechanical interference it can in 
principle occur wherever phase interference phenomena 
exist and are governed by Hermitian eigenvalue prob- 
lems, as in the case of classical electromagnetic waves in 
loss-free media. 

This paper is organized as follows. In section II, we 
present the basic formalism of the Maxwell normal mode 
problem in periodic, loss-free media, discuss the Berry 
curvature of the photonic bandstructure problem, and 
consider the effects of broken time-reversal symmetry. 
In section III, we provide explicit numerical examples of 
bandstructures containing non-trivial topological prop- 
erties, and show the occurrence of edge states along do- 
main wall configurations. Motivated by the numerical re- 
sults, in section IV, we derive a simple Dirac Hamiltonian 
from the Maxwell equations using symmetry arguments 
as the guiding principle, and we show how under certain 
conditions the zero modes of this Dirac Hamiltonian ex- 
hibit anomalous currents along a single direction due to 
the breaking of time-reversal symmetry. It is these zero 
modes that play the role of the "gapless" edge excita- 
tions, as we shall consider in detail. Section V contains 
semiclassical analysis, and we end with a discussion and 
point out possible future directions in section VI. 



II. BERRY CURVATURE IN THE MAXWELL 
NORMAL-MODE PROBLEM 

In this section, we discuss the formal basis of Berry cur- 
vature in the photon band problem. We begin with the 
basic formulation of the photonic bandstructure prob- 
lem, which is somewhat more complicated than the elec- 
tronic counterpart, due to the frequency response of di- 
electric media. We then briefly review the connection 
between Chern numbers, Berry curvature, and the occur- 
rence of gapless edge modes along the boundary where 
Chern numbers of a given band change. 



A. Basic Formalism 

We will be solving the source-free Maxwell equations 
for propagating electromagnetic wave solutions in lin- 
ear, loss-free media, and will ignore absorption, nonlinear 
photon-photon interactions, and other processes which 
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do not conserve photon number. We also assume that 
the constitutive relations, reflecting the response of the 
media to the electromagnetic waves, are given by local, 
but spatially varying tensors with generalized frequency 
dependence. The Berry phase, and the associated Berry 
curvature, are commonly identified with quantum me- 
chanics, but in fact are more generally associated with 
the adiabatic variation of the complex eigenvectors of a 
Hermitian eigenvalue system as the Hermitian matrix is 
varied. 

In quantum mechanics, this Hermitian eigenvalue 
problem is the time-independent Schrodingcr equation; 
in the photonic context, it is the classical eigenvalue 
equation for the normal modes of the Maxwell equa- 
tions. In order to make the correspondence to the stan- 
dard quantum-mechanical formulation of Berry curva- 
ture clearer, we will use a somewhat unfamiliar Hamil- 
tonian formulation of Maxwell's equations, which is ap- 
propriate for loss- free linear media. However, we should 
emphasize that that our results are in no way dependent 
on the use of such a formalism, and are properties of the 
Maxwell equations, however they are written. 

In such a loss-free, linear medium, the coupling of elec- 
tromagnetic modes having different frequencies can be 
ignored, and the electromagnetic fields and flux densities 
X(r, t), X € {D, B, E, H} will be of the form 

X(r,t)=(x*(r,w)e™ + X(r ) u)e- iut ), (1) 

where the quantities X are in general complex-valued 
functions of position and frequency with the property: 

(x(r,w))* = X(r,-u;). (2) 

The dynamics of these fields are governed by the source- 
free Maxwell equations: 

V x E = iuiB , Vxfl = -iul), (3) 
V • D = , V ■ B = 0. (4) 

Consider a single normal mode A propagating at fre- 
quency lo\. For the moment, ignore any internal polar- 
ization or magnetization modes of the medium, and as- 
sume instantaneous, frequency- independent response of 
the dielectric material. In this limit, the permeability 
and permittivity tensors, defined by the relations 

B a (r,LJ x )=fi ab (r)H b (r,LO X ), (5) 
D a (r,uj x ) =e ab (r)E b (r,iu x ), (6) 

are both positive-definite Hermitian tensors and have 
well-defined, positive definite Hermitian inverses e" 1 ^), 
fi~ b (r). Since we have assumed a linear, loss- free medium 
in which photon number is conserved, it is convenient to 
work with a Hamiltonian formalism: the time-averaged 
energy density of the electromagnetic radiation field is 
given by 

H cm (r) =u e (r) + u m (r), (7) 



where 

u e (r) = ^(^,e- 1 (r)£» A ) (8) 

u m {r)= l -(Bi^r\r)B x ). (9) 

Then, if H cm is the spatial integral of the energy density, 
the fields E and H are given by its functional derivatives 
with respect to the divergence-free flux densities D and 
B: 

6H Cttl = J d 3 rE a 5B a + H a 5B a . (10) 

In the local Hamiltonian formalism, the flux density fields 
D(r) and B(r) are the fundamental degrees of free- 
dom, and they obey the following non-canonical Poisson 
Bracket relations: 

{D a (r), 5V)}pb = e a " c V c <5 3 (r - r'). (11) 

This Poisson bracket generates the Faraday-Maxwell 
equations dJSJ: 

*E = {D(r),H™} PB , ^ = {B(r),H^} PB . (12) 

Note that these equations do not generate the Gauss law 
equations J3J|, but merely ensure that any divergences 
V a D a and V a B a are constants of the motion; the Gauss 
laws are additional constraints that are compatible with 
the Faraday-Maxwell equations of motion. 

If internal polarization and magnetization modes of the 
medium are ignored, a discretized form of the electro- 
magnetic Hamiltonian is formally identical in structure 
to that of a collection of real oscillator variables Xi with 
non-canonical Poisson brackets 

{Xi,Xj} P B = Sij, (13) 

where Sij is a real antisymmetric matrix, and the Hamil- 
tonian energy function is 

ij 

where Bij is a real-symmetric positive-definite matrix. It 
is useful to introduce the imaginary antisymmetric Her- 
mitian matrix Ay = iSij] The canonical normal modes 
are given by 

q x ±ip x = (j x )~ 1 ^2uf x x u (15) 

i 

where 7a is an arbitrary scale factor, and where (u^)* 
= u7° , a = ±, which obeys the generalized Hermitian 
eigenvalue equation 

Y, A v u % = ± ^H B iM» ( 16 ) 
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with uj\ > 0, and the orthogonality condition 



XX' 



(17) 



Because of the antisymmetric Hermitian property of 
the matrix Ay, and the positive-definite real-symmetric 
property of the matrix j3y, this eigenproblem has real 
eigenvalues that either come in pairs, ±oj>, or vanish; 
these equations provide a straightforward transformation 
to canonical form only if the generalized eigenvalue prob- 
lem has no zero-frequency eigenvalues, which is only the 
case if Ay is non-singular. 

The coefficients Uix are the analogs of the electromag- 
netic fields E(r,uj) and H(r,oj). It is also useful to in- 
troduce the conjugate quantities 



J2 B 



X<xY 



Saa'SxX' 



(18) 



these are the analogs of the flux densities D(r,oj) and 
B(r,u>), and By encodes the "constitutive relations" be- 
tween analogs of fluxes and fields. 

The Hamiltonian formulation of the Maxwell equations 
indeed presents the difficulty of having a null space of 
zero-frequency eigenvalues: by themselves, the Faraday- 
Maxwell equations have static (zero-frequency) solutions 
B(r) = V/(r), D(r) = Vg(r); the role of the additional 
Gauss law constraints is precisely to eliminate these zero 
modes. The zero-mode problem in the Hamiltonian for- 
mulation is the counterpart of the gauge ambiguity of 
the solutions of Maxwell's equations in the Lagrangian 
formulation. 

In the Maxwell equations, -By becomes the following 
positive-definite 6x6 Hermitian matrix, 



Bi 



»a b \r) 



(19) 



More precisely, this a 6 x 6 block of an infinite- 
dimensional "matrix" that is block-diagonal in terms of 
the spatial coordinate r. (The "A" and U B" matrix no- 
tation is common in the context of generalized Hermitian 
eigenvalue problems, where the positive-definite charac- 
ter of the U B" matrix guarantees reality of the eigenval- 
ues; hopefully the context should distinguish our use of 
the symbol B for such a matrix from the symbol B(r) 
used for the magnetic flux density.) In this continuum 
limit, the antisymmetric Hermitian matrix Ay becomes 
a 6 x 6 matrix block of differential operators: 



A ac = 



ie abc \7 b 
e abc \7 b 



(20) 



This A-matrix can be also be elegantly expressed using 
the 3x3 spin-1 matrix representations of the angular 
momentum algebra, (L fc ) 



abc . 



A = 



it 

L a V a 
-L a V a 



(21) 



From the antisymmetry of A, it again follows that its 
eigenvalues come either in ± pairs, or are zero modes, 
corresponding to static field configurations. Due to the 
presence of a huge band of zero modes (one third of the 
spectrum), the A matrix cannot be written in canonical 
form. 

Using the Poisson-Brackets, we see that the equation of 
motion of the electric and magnetic fields is a generalized 
Hermitian eigenvalue problem of the form 



L a V a 
-L a V a 



E x 



e(r) 
M (r) 



In this formalism, the energy-density of the normal mode 
(time-averaged over the periodic cycle) is 



1 



u(r) = -[E* x HI ] 



e(r) 








E x 
Hx 



(22) 



and the period-averaged averaged energy-current density 
(Poynting vector) is 



f{r)= l -[El HI 




iL a 



-iL a 




H x 



(23) 



For practical real-space-based calculations of the pho- 
tonic normal mode spectrum with inhomogeneous local 
constitutive relations, it is very convenient to discretize 
the continuum Maxwell equations on a lattice (or net- 
work) in a way that they in fact have the matrix form 
(|16f) . where the matrix Ay reproduces the zero- mode 
(null space) structure of the continuum equations, and 
Hij represents the local constitutive relations at network 
nodes (which come in dual types, electric and magnetic). 
In such a scheme, divergence-free electric and magnetic 
fluxes flow along the links of the interpenetrating dual 
electric and magnetic networks, while electromagnetic 
energy flows between electric and magnetic nodes, along 
the links between nodes of the network, satisfying a local 
continuity equation (see Appendix iBll. However, there is 
one further conceptual ingredient that needs to be added 
to the formalism before we can discuss the Maxwell nor- 
mal modes in "non-reciprocal" media which have broken 
time-reversal symmetry. 



B. Frequency dependence of the dielectric media 

The formalism discussed so far treats the constitutive 
relations as static. In general, although we will treat 
them as spatially local, we cannot also treat them as in- 
stantaneous, and must in principle treat the local permit- 
tivity and permeability tensors as frequency-dependent, 
e(r) — > e(r,u), /x(r) — > fi(r,ui). This is because a non- 
dissipative time-reversal-symmetry breaking component 
of these tensors is both imaginary and antisymmetric (as 
opposed to real symmetric) and is an odd function of 
frequency, making frequency-dependence inescapable in 
principle. 
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These effects can on the one hand be treated in a 
Hamiltonian formulation by adding extra local harmonic 
oscillator modes representing local polarization and mag- 
netization degrees of freedom of the medium that couple 
to the electromagnetic fields. The full description of this 
is again a set of harmonic oscillator degrees of freedom 
described by equations of the form (|16|l . On the other 
hand, with the assumption that we are treating the elec- 
tromagnetic modes in a frequency range that is not res- 
onant with any internal modes of the medium (i.e., in 
a frequency range where the loss-free condition is sat- 
isfied), we can eliminate the internal modes to yield a 
purely-electromagnetic description, but with frequency 
dependent constitutive relations. 

The details are given in Appendix El but the result 
can be simply stated. If all oscillator degrees of free- 
dom are explicitly described, the eigenvalue problem for 
the normal modes has the structure (|lfci|> . where B^ 
is positive-definite and real-symmetric. This guarantees 
that the eigenvalues lo\ are real. However, the normal 
modes in some positive frequency range ui\ < u < ui 2 can 
be treated by eliminating modes with natural frequen- 
cies outside that range, to give an matrix-eigenvalue-like 
problem of much smaller rank of the form 

£ A ijU % = ±u x B^{u x )u%, (24) 

3 3 

where Bij (u>) is now a frequency-dependent matrix with 
a Kramers-Kronig structure. The matrix Bij(w) is no 
longer in general real-symmetric, but provided the elim- 
inated modes are not resonant in the specified frequency 
range it instead is generically complex Hermitian. The 
"eigenvalue equation" is now a self-consistent equation: 

Aii«fxb>) = ±«>aM E B~\u>)u%(lj). (25) 

3 3 

This must be solved by varying u> till it coincides with 
an eigenvalue. Unfortunately, while B~^{uj) is Hermi- 
tian and non-singular in the dissipationless frequency 
range u>i < u < u>2, it is not necessarily positive defi- 
nite, so a priori, the eigenvalues lu\(lu) are not guaran- 
teed to be real, except for the fact that we know that 
these equations were derived from a standard frequency- 
independent eigenvalue problem which does have real 
eigenvalues. As shown in Appendix ^ the Kramer- 
Kronig structure of B^ (uj) reflects the stability of the un- 
derlying full oscillator system, giving the condition that 
a modified matrix 

^V)) = |;M^V)) (26) 

is positive-definite Hermitian in the specified frequency 
range, which is sufficient to guarantee reality of the eigen- 
values in that range. Furthermore, the quadratic expres- 
sion for the energy of a normal mode solution is given 
in terms of BZ ((jJ\) rather than B^ (ui\), reflecting the 
fact that the total energy of the mode resides in both the 



explicitly-retained degrees of freedom (the "electromag- 
netic fields" ) and the those which have been "integrated 
out" (the non-resonant polarization and magnetization 
modes of the medium): 

Xi(t) = ^(^u+e'^+cc, cj a >0, (27) 
H = 5£^W«&)*«+. ^=0- (28) 

If the frequency dependence of the positive-definite 
Hermitian matrix Bij (to) is negligible in the range u)\ < 
lu < uj 2 , so Bjj(cj) » B, 3 (w ), with lu = (u>i + w 2 )/2, 
one can replace B^ (u) in (|24H by the positive-definite 
frequency-independent Hermitian matrix B^ (u>o). This 
in turn allows the eigenvalue problem to be transformed 
into the standard Hermitian eigenvalue problem 

H ijW ^ =ljxwI X) , (y>W,wM)=8xx,, (29) 

with scalar product 

(«y)sj;^, (30) 

3 

valid for positive lu\ in the frequency range where B^ (oj) 
»B«(wo), with 

Hij = (B^MAB 1 ^)) , (31) 

and 

3 

This allows well-known Berry-curvature formulas from 
the standard Hermitian eigenproblemA to be quickly 
translated into the generalized problem. It turns out 
that when the full problem with frequency-dependent 
constitutive relations is treated, the standard formula for 
the Berry connection remains correct with the simple re- 
placement Bij(ujo) — > Bij(ujx) ( the Berry curvature and 
Berry phase can both be expressed in terms of this Berry 
connection) . 

C. Berry curvature in Hermitian eigenproblems 

Let Hij (g) be a family of complex Hermitian matri- 
ces defined on a manifold parameterized by a set g of 
independent coordinates g^, \i = 1,...D It is assumed 
that the matrix is generic, so its eigenvalues are all dis- 
tinct; as it is well known, three independent parameters 
must be "fine-tuned" to produce a "accidental degener- 
acy" between a pair of eigenvalues. Thus if the paramet- 
ric variation of the Hermitian matrix is confined to a two- 
parameter submanifold, each eigenvalue 0J\(g) will gener- 
ically remain distinct. Under these circumstances, the 
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corresponding eigenvector is fully defined by the eigen- 
value equation and normalization condition, up to mul- 
tiplication by a unimodular phase factor, that can vary 
on the manifold: 



w\i(g) 



(33) 



This is the well-known U U(1) gauge ambiguity" of the 
complex Hermitian eigenproblcm. Associated with each 
eigenvector is a gauge-field (vector potential in the pa- 
rameter space), called the "Berry connection": 

4 A )( S ) = ^(w A (g)),a A1 w A (g)), = A (34) 

This field on the manifold is gauge-dependent, like the 
electromagnetic vector potential A(r), but its curvature 
(the "Berry curvature"), the analog of the magnetic flux 
density B{r) — V x A(r), is gauge invariant and given 
by 



T${g) = d,A^{g)-d v A^\c 



(35) 



The Berry phase associated with a closed path Y in pa- 
rameter space is given (modulo 2it) by 



exp (-i0 (A) (r)) = exp (-ij^A^g 



W 



(36) 



Ignoring frequency-dependence, the oscillator system 
has 



Hijg) = B 1 / 2 (g)AB 1 ^(g), 



(37) 



where the positive-definite Hermitian matrix B{g) can 
continuously vary as a function of some parameters g, 
but A is invariant. Then converting to the oscillator 
variables gives 

f(uW,B-i(g,u; x )d,u(^ 
4 = Im ' ^7 = 77T" I • (38) 

Here B^ 1 (g) has been replaced by B _1 (<7, lj>,) which 
is the correct result when frequency dependence of 
B~ 1 ig,uj) is taken into account (see Appendix lA"|l . 



D. Photonic bands and Berry curvature 

In the case of periodic systems, the normal modes have 
discrete translational symmetry classified by a Bloch vec- 
tor k defined in the Brillouin zone, i.e., defined modulo 
a reciprocal vector G. For fixed k, the spectrum of nor- 
mal mode frequencies to n {k) is discrete, labeled by band 
indices n, and, as emphasized by Sundaram and Niui in 
the electronic context, the Bloch vector of a wavepacket 
plays the role of the control-parameter vector g. 

In order to compute the Berry curvature of the photon 
band Bloch states, we shall find it convenient to work 



in a fixed Hilbert space for all Bloch vectors k, and we 
do this by performing a unitary transformation on the 
A "matrix" (which becomes the 6x6 matrix of differen- 
tial operators (|2U|I in the continuum formulation of the 
Maxwell equations) as 



A{k, V) = e- lk r A{V)e 



Ak-r 



A(V + ik) 



(39) 



Note that parametric dependence on the Bloch vector k is 
a little different from parametric dependence on param- 
eters g that control the Hamiltonian. After projection 
into a subspace of fixed k, the "^4" matrix also becomes 
parameter-dependent, while (if the constitutive relations 
are taken to be completely local) the "£?" matrix in the 
photonics case is only implicitly fe-dependent through its 
self-consistent dependence on the frequency eigenvalue. 
(Parameter-dependence of the "A" matrix does not af- 
fect the expression (|38fl for the Berry connection.) 

The discrete eigenvalue spectrum uo n {k) is then ob- 
tained by solving the self-consistent matrix-differential- 
cquation eigenvalue problem: 

A(k,V)u n (k,r) =w n (fc)B- 1 (r,w„(fc))M„(fc,r), (40) 

where B _1 (r,tj) is the 6x6 block-diagonal permittivity- 
permeablity tensor diag(e(r, uS), £t(r, u>)), and 
u n {k, r) expifc • r represents the 6-component com- 
plex vector (E n (k,r), H n (k,r)) of electromagnetic 
fields of the normal mode with Bloch vector k and fre- 
quency u) n (k). The eigenfunction satisfies the periodic 
boundary condition u n (k,r + R) = it„(fc,r), where 
R is any lattice vector of the photonic crystal, where 
B-^r + R,u) = B-V.w). 

The transcription of Eq.JJSJ) to the case of periodic 
media then gives the three-component Berry connection 
(vector potential) in fe-space as 



^ n) (fc)=Im. 



M„(fc), J B- 1 K(fc))V^M„(fc) 



u n (k), j B- 1 (w„(A;))m„(A;) 



(41) 

The scalar products in Ea. (|41|l are defined by the trace 
over the six components of u n (k,r), plus integration of 
the spatial coordinate r over a unit cell of the photonic 
crystal. By construction, if a "Berry gauge transforma- 
tion" u n (fc, r) — > u n {k,r)expixn{k) is made, A? n Jk) 
^A? n) (k) + Vl Xn ik). 

In three-dimensional fc-space, the antisymmetric Berry 
curvature tensor ^(fe) = V&AL^fc) - V|-A? n) (fc) can 
also be represented as the three-component "Berry flux 
density" fi!™^(fc) = e a bc^ b k A c ^{k) (the fc-space curl of 
the Berry connection), to emphasize the duality between 
r-space and fe-space, and the analogy between Berry flux 
in fc-space and magnetic flux in r-space, 

If a wavepacket travels adiabatically (without inter- 
band transitions) through a region with slow spatial vari- 
ation of the properties of the medium, so the photonic 
normal-mode eigenvalue spectrum can be represented as 
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a position-dependent dispersion relation u> n (k,r), the 
wavepacket must be accelerated as its mean Bloch vector 
k slowly changes to keep its frequency constant. When 
translated into the language of photonics, the scmiclassi- 
cal electronic equations of motion then become the equa- 
tions of ray optics: 



. a dka 
dt 
dr a 
~~dt 



-h a V a to n (k,r), 
V a k u> n (k,r)+^ b (k,r 



dki, 
~dt' 



(42) 
(43) 



where n oc dr/dt is parallel to the ray path, V a = d/dr a 
and V£ = d/dk a (it is useful to use covariant and con- 
travariant indices to distinguish components of spatial 
coordinates r a from the dual Bloch vector components 
k a ). The Bloch-space Berry curvature f^ b (k,r) controls 
the additional "anomalous velocity^" correction in 143|) 
to the familiar group velocity of a wave packet v®(k) — 
VJJw n (fe), which is active only when the wavepacket is 
being accelerated by the inhomogeneity of the medium. 

Before we conclude our general discussion on Berry 
curvature in photon band systems, we must state the 
constraints inversion and time-reversal symmetries place 
on the tensor T^ h {k). In what follows, we will use the 
Bloch state w n defined in Eg . 1^3211 . If inversion symmetry 
(I) is present, the periodic part of the Bloch state w n (k 
has the following property: w n (k) = w n (—k) whereas 
if time-reversal symmetry (T) is present, w n (k) = 
«>*(— fe). If only (I) is present, it then follows that 
Ff{k) = F% b (-k), whereas if only (T) is present, 
J-"% = — (— fc). If both symmetries are present, then 
the Berry Curvature is identically zero everywhere except 
at isolated points of "accidental degeneracy" , where it is 
not well-defined. When Tf b is non-zero, the phases of 
the Bloch vectors cannot all be chosen to be real. These 
properties will be crucial when we consider the effects of 
various symmetry breaking perturbations on the photon 
band structure. 



E. Topological structure of the photon bands 

The main consequence of having bands with non-zero 
Berry curvature field is that if the path C is closed and 
encloses an entire Brillouin zone, the single- valuedness of 
the state w n requires that 

exp (j> A^dkA = exp (J J dk a A dk b T^j = 1 

or, 



Sc 



dk a Adk b ^ b = 2ttCW, 



(44) 



where Cn^ is an integer, known as the Chern invariant 
associated with the nth band, and have well-known con- 
sequences in the quantum Hall effect: in the integer quan- 
tum Hall effect, where the interactions among electrons 



are weak, the Hall conductance is expressed in terms of 
the sum of all Chern invariants of bands below the Fermi 
levelii: 



(45) 



The gauge structure of the photon band problem out- 
lined above is formally analogous to the local U(l) gauge 
invariance of ordinary electromagnetism. Note that the 
gauge invariance refers to the phase of the of the six- 
component electromagnetic fields as a whole; adding ar- 
bitrary phase fc-dependent phase factors to each field sep- 
arately will in general not preserve the Maxwell equation 
constraint. 

A phase convention can be specified, for instance, 
by arbitrarily choosing real and imaginary axes of the 
phases; the local gauge-dependent phase fields of the elec- 
tromagnetic Bloch states are then represented as two- 
component rotor variables at each point of the Brillouin 
zone. In addition, a gauge choice may be made sepa- 
rately for each band so long as the spectrum remains 
non-degenerate. 

By representing the phase covering on the Bloch mani- 
fold this way, the possibility of the occurrence of topolog- 
ical defects of the gauge field become transparent. Local 
gauge transformations correspond to local smooth defor- 
mations of the rotor variables, and the Chern invariant 
corresponds to the total winding number of theses rotor 
variables along a closed path enclosing the entire Bril- 
louin zone. 

In the case of two dimensional Bloch bands, the de- 
fects of the phase field are point singularities, having a 
zero-dimensional "core" region where a phase convention 
is not well defined, due to quasi-degeneracies with neigh- 
boring bands. It is clear that Bands can have non-zero 
Chern numbers only if time-reversal symmetry is broken. 
Otherwise, the Berry curvature will be an odd function 
of k, and it's integral over the entire 2D Brillouin zone 
vanishes. 

In three dimensions, the defects of the phase field 
are line defects or "vortices" and their stability requires 
quasi- degeneracies to occur along isolated lines in recip- 
rocal space. 

In the photonic system of interest, even if photon bands 
can have non-zero Chern numbers, there can be no Hall 
conductance as given above due to their Bose statistics 
(and hence, to their finite compressibility). However, the 
connection between edge modes and Chern invariants is 
independent of statistics: if the Chern number of a band 
changes at an interface, the net number of unidirection- 
ally moving modes localized at the interface is given by 
the difference of the Chern numbers of the band at the 
interface. We shall consider the problem of how Chern 
numbers can change across an interface in the next sec- 
tion. 

Since the Chern invariant of a band is a topological 
number, it therefore cannot vary smoothly as we vary 
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(b) 

FIG. 1: A representation of the phase fields of the photonic 
Bloch states in a two dimensional Brillouin zone using two 
component rotors. The entire set of six electric and mag- 
netic fields is associated with a single phase at each point in 
the Brillouin zone. The Chern invariant simply represents 
the winding number of this phase along the Brillouin zone 
boundary and is also given by the integral of the Berry cur- 
vature T xy over the two dimensional Brillouin zone. The 
phases in (a) correspond to bands with both inversion and 
time-reversal symmetries, and the phases of the band can be 
chosen to be real every where in the Brillouin zone. For bands 
having non-zero Chern invariants (b), the phase around the 
zone boundary winds by an integer multiple of 2tt and there 
is a phase vortex-like singularity somewhere in the Brillouin 
zone where the Berry connection cannot be defined, due to 
the occurrence of quasi-degeneracies. 




FIG. 2: The number of forward minus the number of back- 
ward moving edgemodes equals the difference of the Chern 
number of the band across the interface. 



some parameter of the periodic eigenproblem. So long 
as a band remains non-degenerate, it's Chern number 
cannot vary. However, if we tune some parameter A of 
the Hamiltonian to a critical value A c such that two bands 
having non-zero Chern invariants touch at some isolated 
point in the Brillouin zone when A = A c , the two bands 
can exchange their Chern numbers at these degenerate 
points; if we tuned A beyond its critical value, the bands 
would emerge with different Chern invariants. Since the 
total Berry "magnetic flux" of all bands remains fixed 



always, if only two isolated bands exchange their Chern 
numbers at points of degeneracy, the sum of their Chern 
numbers must remain invariant. 

Generically, 2D bands with both time-reversal and in- 
version symmetry touch at isolated points of accidental 
degeneracy in a linear conical fashion, forming "Dirac 
cones" in the vicinity of which the spectrum is deter- 
mined by a massless Dirac Hamiltonian: 



H = uj-uj d = v D (Sk 1 a 1 +Sk 2 a 2 ) , (46) 

where vd, is a parameter that gives the slope of the cone 
close to the accidental degeneracy. 




FIG. 3: As we tune some parameter A of the Hamiltonian 
across a critical point where "accidental degeneracies" occur, 
and two bands touch in a linear fashion forming a "Dirac 
point", Chern numbers of bands may be exchanged. 



III. BROKEN T AND I IN PHOTONICS 

In this section, we shall discuss our strategy for con- 
structing photon bands with non-zero Chern invariants, 
and "crural" edge states, whose existence is confirmed in 
the following sections. 

To break time- reversal symmetry in photonics, we shall 
need magneto-optic materials (i.e. a Faraday rotation ef- 
fect). Such materials are characterized by their ability to 
rotate the plane of polarization of light, when placed in a 
magnetic field, and are used in conventional optical iso- 
lators. The amount of rotation per length is known as a 
Verdet coefficient, which depends on temperature as well 
as on the wavelength of light. Materials known to have 
large Verdet coefficients (~ 10°mm~ 1 at wavelengths of 
the order of microns) are the iron garnet crystals such as 
Ho3Fe50i2^i. Due to the breaking of time-reversal sym- 
metry in this materials, the eigenfrequency degeneracy is 
lifted for light characterized by different states of circular 
polarization. 

While such magneto-optic devices employ magnetic 
fields in the direction of travel of the light beam, we shall 
be interested in two dimensional photonic crystals with 
the magnetic fields placed perpendicular to the plane of 
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FIG. 4: In the conventional Faraday effect used in optical iso- 
lators, light travels in the same direction as the applied field, 
resulting in the rotation of its polarization plane. However, in 
the photonic analog of a 2DEG heterojunction, light travels 
in F 



propagation of light, as shown in Fig. We shall call the 
axis perpendicular to the 2D photon bands the Faraday 
axis, and the setup here is reminiscent of a 2D electron 
gas placed in a perpendicular magnetic field. 

Although we now have a means of introducing bro- 
ken time reversal symmetry, we still need a strategy for 
the nucleation of equal and opposite pairs of Chern in- 
variants on bands near points of accidental degeneracy. 
To do this, we choose hexagonal lattice geometry. The 
threefold rotation symmetry of such a system guarantees 
the existence of Dirac points in the Brillouin zone cor- 
ners when both inversion and time-reversal symmetry are 
present; in this case the only irreducible representations 
of the space group of three-fold rotations correspond to 
non-degenerate singlets and degenerate doublets. As a 
simple example consider the case of free photon "bands" 
with dispersion ui = c\k\ in the first Brillouin zone of a 
triangular lattice. The eigcnfrcqucncies of the photons 
are six-fold degenerate at the zone corners. Adding a 
weak periodic perturbation in the constitutive relations 
will lift the degeneracy and the bands will now be either 
non-degenerate or will form degenerate doublets, as de- 
manded by the symmetry of the perturbation. Due to 
the 6-fold rotation symmetry, the doublets are allowed 
to have a linear dispersion close to the zone corners and 
shall be our Dirac points of interest, whereas the non- 
degenerate singlet bands disperse quadratically. We shall 
provide explicit examples of hexagonal photonic band- 
structures having Dirac points in section V. 

While the existence of such Dirac points are virtually 
guaranteed in triangular lattice systems, their stability 
has little to do with lattice geometry. Such points are 
stable in two dimensions only because of the presence of 
time-reversal symmetry and inversion symmetry, when 
the eigenvalue problem is a real-symmetric one: in this 
case it is possible to find "accidental" degeneracies by 
varying just two parameters, according to the Wigner- 
von Neumann Theorem. Thus, if the perfect hexago- 
nal geometry of the constitutive relations is slightly dis- 
torted, the Dirac points will simply move elsewhere in the 
two dimensional Brillouin zone. Provided that such dis- 
tortions are not too strong such that an axis of two-fold 



rotations is introduced, in which case the linear disper- 
sion characteristic of a Dirac point is no longer allowed, 
or if the distortion is so great that the Dirac points meet 
and annihilate at a point of inversion symmetry, Dirac 
points will still exist in the system. 

If, however, inversion or time-reversal symmetry is bro- 
ken in the system, the eigenvalue problem becomes com- 
plex Hermitian, and according to the Wigner-von Neu- 
mann theorem, three parameters are required to ensure 
stability of the Dirac points - in this case, the Dirac point 
degeneracy of the 2D bandstructure is immediately lifted. 
In both cases, the two bands which split apart each ac- 
quire a non-zero Berry curvature field F xy {k). 

If inversion symmetry alone is broken, T X y{k) is an 
odd function of k, as discussed above. While the bands 
do have interesting semiclassical dynamics due to their 
anomalous velocity, they do not have any interesting 
topological properties since their Chern invariants are 
identically zero. 

On the other hand, if time-reversal symmetry alone is 
broken, via the Faraday coupling, the Berry curvature 
field will be an even function of k, and each band which 
split apart due to the Faraday coupling will have equal 
and opposite non-zero Chern invariants. 

Finally, if we can slowly tune the Faraday coupling in 
space, from a positive value, across the critical value of 
zero, where the local bandstructure problem would per- 
mit Dirac spectra, to a negative value, we would gen- 
erate a system of photonic bands with non-zero Chern 
numbers, that get exchanged at the region of space cor- 
responding to the critical zero Faraday coupling. It then 
follows, that modes with exact correspondence to the in- 
teger quantum Hall edge states would arise in such a 
system. In the following section, we shall display this 
explicitly using an example Bandstructure. 



IV. EXPLICIT REALIZATION OF EDGE 
MODES 

An example of a photonic bandstructure with the de- 
sired properties is shown in Fig. |3J It consists of a trian- 
gular lattice of dielectric rods (e a = 14) placed in a back- 
ground of air (e = 1) with a area filling ratio of / = 0.431. 
The authors of ref. I22I in a quest for optimal photonic 
bandgap materials, first studied this system. They com- 
puted the TE mode spectrum and found a full band gap 
in the TE spectrum. We have reproduced numerically 
their calculation and have also computed the spectrum 
for the TM modes. 

The key feature of this particular system which is of 
importance to us are the presence of a pair of Dirac points 
in the spectrum of the TE modes which are well isolated 
from both the remaining TE and TM modes. Each of 
the six zone corners contains the Dirac cone spectrum, 
but there are only two distinct Dirac points, the others 
being related by reciprocal lattice translations of these 
points. In this particular system, the two Dirac points 
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FIG. 5: Photon bands in the k z = plane of a 2D hexagonal 
lattice of cylindrical dielectric rods. Electromagnetic waves 
are propagating only in the x — y plane (Brillouin zone shown 
in the lower right). As in Ref. 1221 the rods have a filling 
fraction / = 0.431, e = 14, and the background has e = 1. 
The band structure contains a pair of Dirac points at the zone 
corners (J). 



FIG. 6: Phase diagram of the photonic system as a function 
of inversion and time-reversal symmetry breaking. In regions 
I and III, the gap openings of both Dirac points are primarily 
due to inversion symmetry breaking, whereas in regions II and 
IV, the breaking of time-reversal symmetry lifts the degener- 
acy of the bands which formed the Dirac point. In all four 
regions, the two bands of interest have non-zero Berry curva- 
ture, but only in regions II and IV do they contain non-zero 
Chern numbers. 



are related by inversion in k-space. 

As we have discussed, a gap immediately opens when 
either inversion or time-reversal symmetries are broken 
in this system. We break inversion symmetry in the sim- 
plest possible way by introducing a slight imbalance in 
the value of the dielectric tensor in the rods at opposite 
ends of the unit cell, and we parameterize the inversion 
breaking by defining the quantity 



Mi = log 



(47) 



where e+ (e_ ) is the value of the permittivity inside the 
rods in the upper (lower) half of the unit cell depicted in 
Fig. El 

To break time-reversal symmetry, we add a faraday 
effect term in the region outside the rods. This is done by 
giving the dielectric tensor a slight imaginary component 
without varying the constitutive relations inside the rods: 



outside rods: 



L (x) 



-iA e- 1 



inside rods: e i7 - 1 (x) = ( e ^ ^ 



We also define a parameter 



M T = A, 



(48) 



(49) 



(50) 



to define the strength of the time-reversal symmetry 
breaking perturbation. 

We first determined the phase diagram of the system 
in the (mj,mT) plane by breaking both inversion and 
time-reversal symmetry, and locating special values of the 
symmetry breaking parameters that result in the closing 
of the bandgap at one or more Dirac points (Fig. [BJ. The 
phase diagram separates regions characterized by bands 
just below the band gap having a non-zero Chern number 
from regions with all bands having zero Chern numbers. 
The boundary between these phases are where the gap at 
one or more of the Dirac points vanishes, as shown in Fig. 
El Since there are two Dirac points, each phase bound- 
ary corresponds to the locus of parameters for which the 
gap at one of the Dirac points closes. Thus, both Dirac 
points close only when both lines intersect, namely at the 
point (raj = 0, my = 0), where both inversion and time- 
reversal symmetries are simultaneously present. When 
inversion symmetry alone is broken, the Berry curvature 
field of Dirac point 1 is equal in magnitude and oppo- 
site in sign of the Berry curvature at the second Dirac 
point. When time-reversal symmetry is broken, on the 
other hand, each Dirac point has an identical (both in 
magnitude and sign) Berry curvature field. In this case, 
the photon bands which split apart at the Dirac point 
each have Non-zero Chern number, which depends only 
on the direction of the Faraday axis (±£). 

We have also studied numerically the frequency gap 
as a function of the time-reversal breaking perturbation 
above and found that so long as the dielectric tensor re- 
mains positive-definite, the gap increases linearly with 
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FIG. 7: The bandgap opened by time-reversal breaking as a 
function of the strength of the Faraday coupling shows that 
the gap is linearly proportional to A. 

e xy (J7J. This will be important when we consider ef- 
fective Dirac Hamiltonians for this problem: as we shall 
see, the fact that exactly at the zone corner, the gap 
rises linearly with Mt is consistent with the spectrum of 
a massive Dirac Hamiltonian with mass Mt- Thus, we 
have shown an example of a bandstructure which con- 
tains Dirac points whose gaps can be tuned using time- 
reversal and inversion symmetry breaking perturbations. 
We can now show the existence of "crural" edge states in 
this system. 

To study edge states in this system, we introduce a 
"domain wall" configuration across which the Faraday 
axis reverses. As we shall now show numerically, (and 
justify analytically in the following section), the edge 
modes that occur along the domain wall are bound states 
that decay exponentially away from the wall while prop- 
agating freely in the direction parallel to the interface. In 
order to study the exponential decay of these modes, we 
glue together N copies of a single hexagonal unit cell 
along a single lattice translation direction R±, which 
shall be the direction perpendicular to the domain wall. 
We treat this composite cell as a unit cell with periodic 
boundary conditions. Since a domain corresponds to a 
certain direction of the Faraday axis, we study a configu- 
ration here in which the axis changes direction abruptly 
across the domain wall from the +z to the — z direction. 

When we consider the spectrum on a torus, there are 
necessarily two domain walls. Furthermore, since many 
unit cells are copied in this system, there are as many 
duplicates of the bands in the enlarged system under con- 
sideration. We study the bandgap precisely at the Dirac 
point as a function of the fractional distance between the 
two domain walls on the torus x (Fig. |SJ for a composite 
unit cell consisting of N = 30 unit cells copied along the 
R± direction. When x = or x = 1, the two domain 
walls are at the same point, and this corresponds to a 
single domain with a single Faraday coupling A. For all 




FIG. 8: In a system with periodic boundary conditions, there 
are necessarily two domain walls separating regions with dif- 
ferent faraday axes. We study the gap of the spectrum at 
the (now non-degenerate) Dirac points as a function of the 
distance x between the two walls. 

other values of x, the "unit cell" comprises a two domain 
system with non-equivalent lengths. In Fig. ^\ the gap 
between the two bands closest to the Dirac frequency de- 
cays exponentially as a function of the distance between 
the two domain walls. We shall show that the exponen- 
tial decay in the gap corresponds to the localization of 
the edge modes along each domain wall. The small gap 
at intermediate values of x, when the two walls are far 
apart corresponds to the fact that each edge mode has 
a small amplitude, and therefore hardly mix with each 
other at those length scales. 
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FIG. 9: The spectral gap between the two bands which split 
apart due to the breaking of time-reversal symmetry. The 
spectrum is computed on a torus for the extended system 
consisting of 30 copies of the hexagonal unit cell. Further- 
more, domain walls, across which the sign of the Faraday axis 
flips are introduced, and the spectrum is plotted as a function 
of the separation x between the walls (see also Fig. |HJ . 

When the domain wall is introduced, translational 
symmetry is still preserved along the direction parallel 
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FIG. 10: The spectrum of the composite system consisting 
30 copies of a single hexagonal unit cell duplicated along a 
direction R±. Both inversion and time-reversal symmetries 
are present, and the Dirac points are clearly visible. While 
the composite system has a spectrum containing many bands, 
only two bands touch at the Dirac point. The dispersion is 
computed in k space along the direction parallel to the wall. 



FIG. 12: Same system as above, but with a domain wall in- 
troduced corresponding to maximum separation of the walls 
on the torus. The two additional modes present in the gap 
correspond to edge modes with a "free photon" linear disper- 
sion along the wall. There are two modes, since across the 
domain wall, the Chern number of the band just below the 
band gap changes by 2. 
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FIG. 11: The same system as above, but with broken time- 
reversal symmetry without a domain wall. There is a single 
Faraday axis in the rods of the entire system. 



to the wall, and the states of the composite system of 
30 unit cells can be labeled by fey, Bloch vectors in the 
direction parallel to the wall. Figures ITTil 1111 and 1 1 21 con- 
sist of a spectral series of a system without any broken 
time-reversal symmetry (Fig. 1100 . with uniformly bro- 
ken time- reversal symmetry fFig lll(l . and a domain wall 
configuration (Fig. 112(1 for the 30 unit cell composite sys- 
tem. The bands are plotted along a trajectory in fe-space 
in the fen direction which contains the two distinct Bril- 
louin zone corners. It is clear that in the Domain wall, 
there are two additional modes formed in the bandgap 
that arose from the Faraday coupling. Since the domain 



walls are duplicated on the torus, the spectrum of edge 
modes will also be doubled; in Fig. ^1 only the two non- 
equivalent modes are shown. Each mode in the band gap 
has a free photon linear dispersion along the direction of 
the wall; moreover, both have positive group velocities, 
and therefore propagate unidirectionally. 

To be certain, however, that these "chiral" modes are 
indeed localized near the interface, we have numerically 
computed (M(r)|B _1 |u(r)), the electromagnetic energy 
density (the B matrix, defined in section II, is not to 
be confused with the magnetic flux density), the pho- 
ton probability density in real space. We have computed 
this quantity along with all the spectra of the composite 
system using the real space bandstructure algorithms de- 
scribed in Appendix iBl As shown in Fig. ^1 the energy 
density is a gaussian function, peaked at the position of 
the domain wall, decaying exponentially away from the 
wall. From this calculation, we extract a localization 
also approximately 5 unit translations in the direction 
perpendicular to the interface. 

We have therefore shown here using explicit numeri- 
cal examples that photonic analogs of the "chiral" edge 
states of the integer quantum Hall effect can exist along 
domain walls of Hexagonal photonic systems with broken 
time-reversal symmetry. We have studied the unphysical 
case in which such domain walls are abrupt changes in the 
axis of the Faraday coupling. However, due to the topo- 
logical nature of these modes, a smoother domain wall 
in which the Faraday axis slowly reverses over a length 
scale much larger than a unit cell dimension would also 
produce such modes. The most important requirement 
for the existence of these modes, is that at some spa- 
tial location, the Faraday coupling is tuned across its 
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critical value. How this particular tuning is effected is 
irrelevant. In the following section, after deriving the ef- 
fective Hamiltonians for these modes, we shall consider 
a smoothly varying Faraday coupling, which corresponds 
to an exactly soluble system, and shall show the evolution 
of these modes as the smoothness of the Faraday coupling 
is varied towards the step function limit considered here. 

0.018 I 1 7v 1 

0.016- / \ 

j? 0.014 - / \ 

w 0.012 - / \ 

I 0.01 - / \ 

>, 0.008 - / \ 

tn \ 

n 0.006 - / \ 

CD / \ 

h 0.004 - / \ 

0.002 - / \ 

o I — ^ 1 

-10 -5 5 10 

Unit cells 

FIG. 13: The real-space electromagnetic energy density pro- 
file associated with the edge modes in Fig. 1121 plotted as a 
function of the direction perpendicular to the domain wall 
(and "integrated" over the direction parallel to the interface) 
and fit to a gaussian profile. The integrated energy density 
depicted here plays the role of the photon probability density 
which confirms that light is confined to the interface. 



V. MODEL HAMILTONIAN APPROACH 

The crucial feature exploited in the previous sections 
was the possibility of tuning bandgaps at Dirac points 
by adding time-reversal breaking perturbations. Before 
adding these perturbations, the linear conical spectrum 
at these Dirac points are governed by two dimensional 
massless Dirac Hamiltonians, and time-reversal or inver- 
sion symmetry breaking perturbations contribute mass 
terms to the Hamiltonian. In this section, we shall 
construct these Dirac Hamiltonians starting from the 
Maxwell equations for two dimensional photonic systems 
with hexagonal symmetry. 

To motivate Dirac Hamiltonians in photonic systems, 
we begin this section by considering a "nearly-free pho- 
ton" approach in which a two dimensional "free photon" 
spectrum consisting of plane waves is perturbed by a 
weak periodic and hexagonal modulation of e(r). Due to 
the underlying symmetry of the perturbation, the plane 
waves mix in a manner to generate Dirac points in the 
zone corners of this system. We then consider the effect 
adding time-reversal and inversion symmetry breaking 
perturbations in this system and derive an expression for 
the Dirac mass. Having motivated the Dirac points, we 



revert to our photon band problem and derive expressions 
for the Dirac mass in these systems. 

In analogy with the "nearly-free electron" approxima- 
tion, we consider the photon propagation problem in the 
weak-coupling regime, in which the dielectric properties 
of the medium act as a weak perturbation. We solve 
the Maxwell normal mode problem for Bloch state solu- 
tions, and work out corrections to the free photon dis- 
persion relations in the Brillouin zone boundaries. We 
shall assume continuous translational invariance in the z- 
direction, and study the propagation of electromagnetic 
waves in the x-y plane. 

The free photon constitutive relations are isotropic and 
uniform in the plane: 

s ° = (o\°.)- 0») 

We consider the "free photon bands" in the first hexag- 
onal Brillouin zone depicted in Fig. [21 Let Gi, i = 1, 2, 3 
be the three equal-length reciprocal lattice vectors each 
rotated 120° with respect to one another. The hexag- 
onal zone corners correspond to the points ±Ki, where 
K x = (G 2 - G 3 ) /3, etc., and \K\ = \G\/y/3. At each of 
the zone corners, the free-photon spectrum is six-fold de- 
generate with too = CqK. In two dimensions, the modes 
decouple into TE (E x ,E y ,H z ), and TM (H x ,H y ,E z ) 
sets, and we shall focus only on the TE modes and con- 
sider the 3-fold TE mode symmetry at the zone corners 
(the TE and TM modes do not mix in 2 dimensions). 
The eigenvalue equation for the free photon plane wave 
modes at the zone corners is A|tto) = ujoBq 1 ]^) , or 
equivalently, B^ 2 AB^ 2 \z Q ) = Wol^o) an d the states 
\zo) = Bq 1/2 \u ) satisfy (z { Q X,> \z ( X } ) = 6\\>. 

A Gi 



K, 




K 3 G, 



FIG. 14: In the weak coupling approach, the free photon TE 
mode plane waves are perturbed by a periodic modulation 
in the permittivity. The plane wave frequency at the three 
equivalent zone corners (Ki, i=l,2,3) is lifted by the permit- 
tivity in "fc • p" perturbation theory into a non-degenerate 
singlet and a degenerate doublet. 

Next, keeping the uniform isotropic permeability fixed, 
we add a weak periodic perturbation to the permittivity 
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of the form 

Asr 1 = ( eoA ^ (r) J), (52) 

with 

3 

V G (r) = 2 cos {G n -r). (53) 

n=l 

After this perturbation is added, the TE and TM 
modes no longer remain degenerate; while the TM modes 
remain 3-fold degenerate at the zone corners at the fre- 
quency ui = cq\K\, the TE modes split apart into a 
singlet and a degenerate doublet. We now determine 
the splitting to leading order in A with within a weak- 
coupling "nearly-free photon" approach. 

With the periodic perturbation, the eigenvalue prob- 
lem is 

A\u)=u;(B 1 + XB^)\u), (54) 

which is equivalent to 

Bl /2 (A - XujB^ 1 ) Bl /2 \z) = (cj + Su) \z). (55) 

The energy splittings are worked out in degenerate per- 
turbation theory (see Appendix IH)l as 

— - = — Hz„\B^ 2 B 1 x B^ 2 \z n ) 

UJQ 

= -A(u„|B 1 ; 1 |m„), 

where \z n ) are appropriate combinations of the three 
free photon plane-plane waves that diagonalizes the pe- 
riodic potential. These states are obtained by requir- 
ing them to be invariant under 3-fold rotations in the 
plane. Instead of writing the fields in the coordinate 
basis, It is convenient to use a redundant basis of the 
three vectors (e iK ^ r , e iK2r , e iK3 ' r ), with J2n K n = 0, 
and Ki ■ Kj = —K 2 /2, i ^ j. In this basis, the magnetic 
field of the TE modes is written as (77 = e 27 ™/ 3 ): 

iff = (1,1,1), (56) 



ff 2 * = (l,T?*,77), (57) 

and 

HI = (1,77,77*). (58) 

The corresponding electric flux densities are easily ob- 
tained: 

d\ = -(z x K u z x K 2 ,z x K 3 ), (59) 



Dl = -(zxK 1 ,r ] *zxK2,T]zxK 3 ), (60) 



D!<= -(zx K llV zx K 2 , V *zx K 3 ) 7 (61) 

and 

!*> = (#!)• ( 62 ) 

Clearly, these are the plane wave solutions that satisfy 
Maxwell equations and transform appropriately under 3- 
fold rotations in the plane. We are therefore led to the 
simple result that the splitting at the zone corners due 
to the mixing of the three plane waves is related to the 
integral over the unit-cell of the electric fields and the 
periodic potential, which is a traceless, real-symmetric 
3x3 problem. It is easy to see that the problem is 
traceless because diagonal terms of the form (ui\B\\ui) 
vanish identically since m are plane waves. 

To leading order in A, the three photon bands 
split to form a singlet band at frequency ujq = 
co|if|(l + A/2 + 0(A 2 )) and a degenerate doublet at fre- 
quency 

lj d = cq\K\ (l-A/4 + 0(A 2 )) . (63) 

Exactly at the zone corners, the singlet and doublet 
states above diagonalize the perturbation in Eq. I|52() . To 
leading order in A and 5k = k — Ki , the deviation in the 
Bloch vector from the zone corners, the states \z2(5k)) 
and \z 3 (5k)), (where \zi(6k)) = exp(i<5fc • r)\zi)), which 
are degenerate at Sk = mix and split apart linearly as 
a function of \Sk\, forming a "Dirac point". To leading 
order, the Dirac point doublet does not mix with the sin- 
glet state \z%(Sk)). The effective Hamiltonian governing 
the spectrum of the doublet, to leading order in Sk is a 
2D massless Dirac equation: 

u±(Sk) = ujd ± v D (5k x <r x + 5k y <T y ) , (64) 

where vd = c /2 + 0(A), and a % are the Pauli matrices 
written in the subspace of the doublet states. The linear 
dispersion of the doublet in the neighborhood of the zone 
corners is immediately obtained by solving Ea. H64|) : 

uj = lue> ± V[)\5k\. (65) 

The singlet band's frequency remains unchanged to 
leading order in 5k: oj {5k) = uj + 0(\5k\ 2 ). Thus, we 
have shown that the periodic modulation of the permit- 
tivity having 3-fold rotational symmetry gives rise to a 
quadratically dispersing singlet band and a "Dirac point" 
with linear dispersion. 
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FIG. 15: Spectrum of photon dispersion in the vicinity of 
the zone corners. We have arbitrarily set A < so that the 
singlet band has a lower frequency than the Doublet. Free 
photon spectra are given by dashed lines. Away from the zone 
corners, the free spectrum is not affected to leading order in 
A. 

Next, we add a Faraday term, with an axis normal to 
the a;y-plane, to the permittivity tensor e xy — —e yx — 
ieorj(r,u}), where 

r)(r,w)=r} (w)+r)i(uj)VG(r). (66) 

Both r)o(Lu) and rjiito) are odd functions of u>. In 
the limit that the Faraday coupling is much weaker in 
strength than the periodic modulation, \r)o\, \r)i\ <C |A| <C 
1, the mixing between the non-degenerate singlet state 
and the doublet remains negligible, and the energy of 
the singlet state is unaffected by the Faraday perturba- 
tion. However, the doublet states split apart at the Dirac 
point. Using the expression for the Dirac point splitting, 
derived in the Appendix [U] we find that the splitting of 
the doublet at the zone corner is given by 

uj± - lo d = ±v d k, k = \K\ Q?7i(wd) - 3Xri (uj D )^ . 

(67) 

Away from the Dirac point (but still close enough to the 
zone corners so that the "nearly-free photon" approxima- 
tion for the three plane wave states remains valid), the 
doublet bands acquire a dispersion 

u = uj D ±v D (\5k\ 2 + k 2 ) 1/2 , (68) 

which is the spectrum of a 2D massive Dirac Hamilto- 
nian: 

uj±(6k) =uj d ±v d {Sk x a x + 5k y <r y + k<t z ) . (69) 

The Dirac points that occur in the "nearly-free pho- 
ton" approximation are not isolated points of degener- 
acy, since, away from the zone corners, the two bands 
which formed the Dirac point merge together to resume 



their original free-photon form. Consequently, the type 
of modes studied in the previous section cannot be re- 
produced using this type of weak-coupling expansion. 

However, we can gain understanding by suppose that 
we have the exact solutions of the electromagnetic Bloch 
states and eigenfrequencies of a system containing iso- 
lated Dirac points, such as the one studied numerically 
in section IV. We can use precisely the same weak Fara- 
day coupling approximation to work out the splitting of 
the Dirac point with a Faraday term. Assuming we are 
given example photonic bandstructures of long hexago- 
nal systems with k z = 0, which contain only isolated 
Dirac points, a weak Faraday coupling would split apart 
the bands that formed the Dirac point, and the split- 
ting is identical to that in (|68J) . Suppose that the two 
bands having a Dirac point, otherwise form a PBG with 
a gap A 3> vdk (as in the case of the numerical exam- 
ple given in the previous section. In this case, since the 
Faraday term removes all points of degeneracy, the now 
non-degenerate bands have a well-defined Berry curva- 
ture field 

T ± (5k) = ±± K (\6k\ 2 + K 2 y 3/ \ (70) 

which decays rapidly away from the Dirac point, and 
contributes a total integrated Berry curvature of ±tt. 
Since there are two non-equivalent Dirac points in the 
hexagonal geometry under consideration, the net Berry 
curvature of the system is the sum of the contributions 
from each Dirac point. If, as in the case under consider- 
ation, spatial inversion symmetry is preserved, but time- 
reversal symmetry is broken, the Berry curvature fields 
at each Dirac point of a given band add, giving total 
Chern numbers ±1 for each of the split bands. However, 
if time-reversal symmetry were preserved, and inversion 
symmetry breaking caused the gap to open, the Berry 
curvature field of each Dirac point for a given band are 
equal in magnitude but opposite in sign, and the Chern 
number would vanish. 

As before, to get unidirectional edge modes of light in 
this system, the Faraday coupling must be tuned across 
its critical value ri(r,u) — 0. To do this, we consider a 
Faraday coupling that varies slowly and adiabatically in 
space, we shall assume negligible frequency dependence 
of the Faraday coupling, and we shall parameterize the 
local value of the Faraday coupling by a smoothly varying 
function k(t), which is positive in some regions and neg- 
ative in other regions of the 2D plane perpendicular to 
the cylindrical axis of the hexagonal array of rods. Due 
to the adiabatic variation of k(t), each point in space is 
characterized by a local bandstructure problem, and the 
splitting at the Dirac point is given again by the expres- 
sion in (|68|) . but with the local value of n. In this limit, 
the smooth variation of n{r) leads to a 2D Dirac Hamil- 
tonian with a adiabatically spatially varying mass gap. 
At all points where n(r) = 0, the local bandstructure in 
the vicinity of the Dirac point is the massless 2D Dirac 
Hamiltonian; provided that \n(r)\ <C A, the PBG, the 
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spectrum far away from the Dirac points is unaffected by 
n(r). In what follows, we assume that when k = 0, our 
bandstructure contains Dirac points which are formed 
by two isolated bands in a PBG region having no other 
points of degeneracy. 

We neglect the mixing between modes at different 
Dirac points, and consider the situation in which k(t) 
vanishes along a single line ( x = for instance), and we 
assume translational invariance along the direction par- 
allel to the interface (y— direction). As before, we con- 
sider the degenerate perturbation problem of the normal 
modes close to the Dirac point. Now, however, the coef- 
ficients of the degenerate solutions of the Maxwell equa- 
tions are spatially varying quantities. Let |w CT (±fe£))), 
a = ±, be the degenerate solutions (i.e. the periodic 
parts of the photon Bloch state wave functions) at a pair 
of Dirac points when k = 0. With the local variation, we 
take spatially varying linear combination of these Bloch 
states 

u(k D , r) = y exp (±ik D ■ r) u a {±k D ,r), (71) 

a,± 

and arrive at the fact that the local value of the splitting 
of the two bands at fen is 

uj+{k D ) -LJ-{k D ) = 2n{r). (72) 

In the neighborhood of the Dirac point, the degener- 
ate perturbation problem gives us a 2D massive Dirac 
Hamiltonian, with 5k x replaced by the operator — i\7 x in 
the position representation, since translation symmetry 
in the x-direction is broken by n(x). We thus obtain an 
expression of the form vdK\^) = 6uj\?p), and 

K = -icr x V x + Sk ll( T y + k(x)(t z . (73) 

The Bloch vector in the y-direction, which remains con- 
served due to the preservation of translation invariance 
along this direction, is ko y + 5k\\. 

For the particular choice of k(x) — k°° tanh(x/£), 
£ > 0, (where k°° is the asymptotic value of the Dirac 
point splitting at distances 3> £ from the interface), the 
problem is exactly solvable, since the Dirac Hamiltonian 
K, when squared, becomes a ID Schrodinger Hamilto- 
nian K 2 corresponding to the integrable Poschl- Teller 
HamiltonianSi. 

To see how this comes about, we explicitly work out the 
operator K 2 , making use of the anti-commuting property 
of the Pauli matrices {a a ,a b } = 2S ab : 

K 2 - 5k 2 = -Vl + k(x) 2 - a v n'. (74) 

The spatially varying Dirac mass term that changed 
sign across the interface becomes a "potential well" with 
bound states given by2i 

Ua{5k\\) = lo d + s K v D 5k\\, s K = sgn(K°°) (75) 

/ \ Y l 2 

u n ± = uj d ±v D (Sk 2 + K 2 ) , n > 0. (76) 



where \n n = 2n|K°°|/£, n < |k°°£/2. In the n — mode, 
light propagates unidirectionally, with velocity , in the 
direction parallel to the wall. All other bound modes 
are bidirectional modes. The numerical example of a 
Dirac mass studied in the previous section that changed 
sign abruptly, as a step function, has the ID Schrodinger 
problem in an attractive delta function potential as its 
square. Consequently, as we have seen, the model per- 
mitted for only a single bound state, corresponding to 
the unidirectional mode. 




-20 -15 -10 -5 5 10 15 20 
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FIG. 16: Spectrum of the integrable Poschl- Teller model. 
With the exception of the zero mode, all bound states corre- 
spond to bi-directionally propagating modes localized at the 
interface where the function n(r) = 0. The zero mode, on the 
other hand, is unbalanced, and furthermore, it corresponds to 
unidirectional propagation. 

For the generic case, the second order differential equa- 
tion for the n > bound states can not be solved ana- 
lytically. However, a formal solution for the zero mode 
eigenfrequency can be obtained, as it is obtained by solv- 
ing a first order equation, as we now discuss. Starting 
from the Dirac equation for the more general case 

v D (-ia x V x + -icr y V y + k(x)(t z ) \tp±) = 5u)\%p±), 

(77) 

by definition, the "zero mode" has the free photon dis- 
persion along the direction parallel to the wall, which 
implies that the function \tjj) oc exp(i5kny). We are thus 
left with the equation 

(-i<r x V x + k{x)<j z ) |^±) = 0. (78) 

Multiplying both sides with cr x , we arrive at the following 
first order differential equation : 

(V a + k(x)0 = 0, (79) 

which has as it's formal solution 

\ip±) = exp \i5k\\y + a J dx' k(x')J \<j>±(a)) , (80) 

where (T v \cj>±(a)) — a\<f)±(a)). Although there are for- 
mally two solutions for the zero mode, corresponding to 
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a = ±1, only one can occur; the other is not normal- 
izable and thus cannot represent a physically observable 
state. 



This make the integral for </>(fc 2 ) trivial (it becomes ex- 
pressible in terms of a simple Hilbert transform), and 
the asymptotic small-fc 2 form H84[l remains valid for all 
values of fc 2 in the domain of the function: 



VI. SEMICLASSICAL ANALYSIS 

Now let the "Dirac mass" term that opens the photonic 
band gap be a slowly varying function k{x) that changes 
monotonically (and analytically) from — fco at x = — oo 
to fco at x — +00. The photonic spectrum of modes 
with wavenumbers k = kjy + 5k near the "Dirac point" 
k[), and which become doubly-degenerate at kp, is an 
adiabatic function of x: 



ui(x,6k x ,5ky) — u> d ± v d ($ky + k(x , Sk x 
k(x,Sk x ) 2 — 5 k 2 + n(x) 2 , 



1V2 



(81) 



where vo > is the "Dirac speed". For k(x,Sk x ) 2 < fc 2 , 
the modes are evanescent as x — > ±00, so are localized 
on the wall. In the x — 5k x plane, the contours of con- 
stant k(x, Sk x ) 2 < fc 2 are simple closed curves, enclosing 
a finite dimensionless area 0(fc 2 ), given by 



</>(fc 2 



= 2 



dx (fc 2 - k(x) 2 ) 



1/2 



(82) 



where x_(fc 2 ) < a;+(fc 2 ) are the two "turning point" so- 
lutions of k(x±) 2 = k 2 . Since k(x) is assumed to be 
monotonic, this can be written as 



4,{k 2 ) = : 



|fe| 



dy (fc 2 - y 2 ) 



2\!/2 



1 



1 



K' + (y 2 ) K'_{y 2 ) 



(83) 



x ± (k 2 ) 



Note that this transformation has turned 4>{k 2 ) into a 
signed area, where sgn(0) = sgn(fco), which is indeed 
the correct form (the function 4>(k 2 ) vanishes as fco — > 
0, when its domain fc 2 < fc 2 , shrinks to zero). In the 
limit fc 2 — * 0, x±(k 2 ) — * xq, the formal location of the 
interface. Then n'j_(k 2 ) — > k'(xq), and ^(fc 2 ) vanishes as 



0(fc 2 ) 



^fc 2 



0. 



re'(xo) 

It is very instructive to examine the special case 
k(x) = fcotanh(a(x — xo)), 
which is integrable. In this case, 

k'(x) — afcosech 2 (a(a; — Xo)), 
k a sech 2 (a(x±(k 2 ) — x )) 



fc 2 



fco 



(84) 

(85) 

(86) 
(87) 



Thus the explicit dependence on x±(k 2 ) can be elimi- 
nated, and 



«±(fc 2 ) = 



fc 2 



(88) 



0(fc 2 ) 



nk^ 
afco 



(89) 



Then the frequency of the interface mode with wavenum- 
ber 5k y = 6k\\ along the interface can be expressed as 



1/2 



xko4>\/TT. 



(90) 



A standard "semiclassical" analysis of interference effects 
on a light ray trapped in a "waveguide" at an interface 
would conclude that the "quantized" values of <f> corre- 
sponding to interface modes were 



b n = 2nn + 7, 



(91) 



where 7 is a "Maslov phase" , usually it. In this case, com- 
parison with the exact solution of the integrable prob- 
lem confirms that this problem instead has a vanishing 
"Maslov phase" 7 = 0. This can be attributed to an un- 
derlying "Z2" Berry phase factor of —1 (Berry phase of 
7r) for orbiting around the degeneracy point at [x— xq, k x ) 
= (0,0). 

We then conclude that the interface modes at a slowly- 
varying interface are in general given (for small #fcii) by 



wo(<5fc||) = ud + «DSgn(ko)£fc| 
w n ±{5k\\ 

<t>(k 2 n 



u>d ± vd (jik 2 + fc^ 



1/2 



,n > 1, 



= 27T71, fc 2 < fcp. 



(92) 



The unidirectional "zero mode" persists however sharp 
the interface is; the bidirectional modes with n > 1 must 
obey 2~Kn < (/'(fco), which has fewer and fewer (and even- 
tually no) solutions as the width of the interface region 
shrinks. In the special case of the integrable model 18511 . 
this spectrum is exact for small 8ku without any condi- 
tion that the wall is slowly varying. 



VII. DISCUSSION 

The occurrence of zero energy modes in the 2D Dirac 
Hamiltonian is well known and represents the sim- 
plest example of a phenomenon known as the "Chiral 
anomaly". The crucial feature, namely, the occurrence 
of interfaces where the Dirac mass gap changes sign, cor- 
responds to tuning our photon band problem across a 
critical point using a Faraday effect. 

We have shown that analogs of quantum Hall effect 
edge modes can exist in photonic crystals whose band 
gaps can be tuned by a Faraday coupling. The crucial 
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new feature we present here here is that photonic sys- 
tems can have bands with non-trivial topological proper- 
ties including non-zero Chern invariants. These in turn 
can be varied in a controlled manner to yield unidirec- 
tional ("crural") edge modes. The edge modes are ro- 
bust against elastic back-scattering since they are states 
which are protected by the underlying 2D band structure 
topology. However, they are not robust against photon 
number non-conserving processes, such as absorption and 
other non-linear effects. We believe that this could be an 
entirely new direction in "photonic band structure engi- 
neering" due to the absence of scattering at bends and 
imperfections in the channel. 




D.O.S D.O.S 

(a) (b) 



A practical realization of such one-way transmission 
channels in photonics will necessarily have to deal with 
the problem of finding a magneto-optic material with a 
strong enough Faraday effect to confine the light close 
to the interface. Furthermore, in a practical design, the 
problem of TE/TM mode mixing when light is confined 
in the direction perpendicular to the 2d system will have 
to be addressed. A practical design could, for instance, 
make use of PBG materials to confine light in the z- 
direction. Although there are many obstacles to the re- 
alization of such interesting effects in photonics, none of 
them arc fundamental, and we believe that these unidi- 
rectional channels could have potentially useful techno- 
logical applications which could in principle be realized 
someday through "bandstructure engineering" . 



FIG. 17: We have shown in this paper that although the 
photon bands (b) cannot be filled as in the electronic case 
(a), they can have no analog of the bulk quantum Hall ef- 
fect. However, the Chern number is a topological invariant 
of Bloch states independent of the constituents. With the 
Faraday term, we are able to tune the system such that the 
total Chern number below a photonic band gap changes across 
in interface, which gives rise to unidirectionally propagating 
edge modes of photons localized at the interface. These modes 
are direct analogues of the "Chiral" edge modes of electronic 
systems which occur at interfaces between two regions hav- 
ing different total Chern invariants below the Fermi level (i.e. 
with different Hall conductances). 
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APPENDIX A: FREQUENCY DEPENDENCE OF 
THE DIELECTRIC MEDIA 

In this section, we shall provide the details of the gen- 
eralization of the normal mode problem to include the 
frequency-dependent response of the media outlined in 
section II. We shall couple the electromagnetic fields to 
harmonic oscillator degrees of freedom of the medium. 
Defining <j>i a and ~Ki a (i = 1, • • • N, a = e, ji) to be a set 
of N independent canonically conjugate oscillator coordi- 
nates and momenta respectively, which represent internal 
polarization and magnetization modes, we consider the 
total Hamiltonian 



H = H* m +Y,H a > = £,/*), (Al) 

a 

where, for instance, 

W = ]T D a (a je (r) Q 7r 4£ (r) + /?? (r)<Mr)) 

i 



The first term above represents the local coupling be- 
tween the electric fluxes and the polarization modes, 
whereas the second term represents the energy of the 
oscillators themselves. A similar equation exists for the 
magnetization degrees of freedom coupled with the mag- 
netic fluxes. The Hamiltonian, as stated in Ea. ljAlfl . is 
real-symmetric and positive-definite, and therefore, its 
eigenvalues are real. The electric and magnetic fields are 
obtained by varying the Hamiltonian with respect to the 
associated flux densities: E a (r) = 5H/SD a (r), H a {r) — 
6H/SB a (r) 



and similarly for the field H a . The time-evolution of the 
oscillator modes are obtained from the Hamilton equa- 
tions of motion (letting dt4>na = —i^>4>ncr, etc) 



iujir ne (r) = 



sn 



Sn ne (r) 

sn 



= u Jne Tr n (r)+^(r)D a (r)(A3) 
w„ e <Mr) + a a ne (r)D a (r).(AA) 



We invert this equation to solve for the oscillator co- 
ordinates and momenta in terms of the fluxes: 



0ne(r) 
7T„e(r) 



LO n ILO 
— lUJ LO n 



D a {r). 
(A5) 

By substituting Ea. (|A5|l into the expression for the 
electric field <|A2|1 . we obtain a correction Se7 t , (r lo) to 
the permittivity tensor coming from the oscillator modes: 



where 



y, f T« 6 (r)(w + u n ) - r*|(r)(w - u n ) 



n 



(A6) 



rf (r) = « e (r) - il3 a ne {v)) (at(r) + i0 b ne (r)) . (A7) 

Finally, the correction term above to the permittivity 
is expressed in Kramers-Kronig form as 



(A8) 



The same formal manipulations occur in the frequency 
dependence of the magnetization modes; in the end, the 
constitutive relations are given by a tensor B(r, lo) de- 
fined by 



B(r,Lj) 



e 1 (r, lo) 







fi- l {r, u>) J ' 
which is written in Kramers-Kronig form as: 

£ — ' \L0 — UJr> LO + LOn 



(A9) 



E a {r)=e-J(r)D b (r) 



^K f (r)i £ (r)+ffi ( (r)i„ £ (r)), 

n 

(A2) 



(A10) 



The first term, S a b(r) = lim a ,_ >00 B(r, lo) is the same 
tensor defining the Hamiltonian in Ea. (|14f> . In the zero 
frequency limit, 

B ab (r, 0) = S ab (r) - £ ( hM+lM ) . (A11) 

Stability of the medium imposes the following impor- 
tant constraint: 



B(r, 0) > 0. 



(A12) 
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Eliminating S a b m Eo. l|A10|) using Eq. i|All|) . we get to be positive-definite in lossless frequency ranges is 



SB{ui) = J2 



u n (w - U) n ) 



u n (u) + uj n ) 



where 8B(uj) = B(uj) — B(0). Whereas B(lu) is not a 
positive-definite matrix, the quantity which is guaranteed 



because 



B(oj) = B(uj) - uj—B(uj) > 0, 

OUJ 



(A13) 



J 



B{uj) = B(0) 



E- 



to — uj,, 



r 



Ul 



UJ + 0J T 



(A14) 



and r„, r*, and B(0) are all positive-definite tensors. 

Although B{oj) is not positive-definite, we will be in- 
terested in cases where 



Deb(B{w)) = 0. 



(A15) 



When this condition is satisfied and B(uj) has no zero 
modes corresponding to metallic conditions, there is a 
well defined inverse tensor B~ 1 (uj) 



B-\r, u) = 



e(r, uj) 







^(r, uj) 



(A16) 



From the stability condition stated for B(uj), there exists 
a similar condition for _B _1 (w): 



B — uj-^—B = B { B- 1 



ouj 



where we have made use of B~ 1 B = 1 and 
d/dw = 0. Supplementing the inequality above 

with the condition in Ea. (|A15|i . we obtain 

B~\u) = ^- (ojB-^lj)) > 0. (A17) 
ouj 

The eigenvalue problem is solved for each value of the 
bloch vector k in the first Brillouin zone, and The formal 
strategy for obtaining the energy eigenvalues is to solve 
A\u n (k)) = uj n {k)B~ 1 (uj{k))\u n {k)) , and then to vary 
uj until it coincides with a frequency of an eigenmodc. 
The stability condition (see Ea. (|A17|l ) guarantees that 
such a prescription enables us to find the entire spec- 
trum in a lossless range of real frequencies, where B^ 1 is 
Hermitian. 

Indeed, if we consider for the moment the Hermitian 
problem 

(A - ujB-^uj)) \u n ) - \ n {uj)\u n ), (A18) 
and vary oj to find the zero modes 

X n (uj) = 0, (A19) 



the stability of such a prescription is guaranteed only if 



duj 



< 0, 



(A20) 



so that the eigenvalues are monotonically decreasing 
functions of uj. But from first order perturbation the- 
ory, we know that the requirement above is satisfied only 
if 



(u n \— (ujB-^uj)) \u n ) > 0, 



(A21) 



which is precisely equivalent to the condition in Eq.( 
fATfjl . 

When we eliminate the internal oscillator (polari- 
ton) modes and explicitly substitute the expressions in 
Ea. (|A5p into the total Hamiltonian, Eq. IjAlfl . we ob- 
tain the following quadratic form that involves only the 
electromagnetic flux densities: 



H 



1 



\ " 6-1 



2 ^ B v 



(uj)(u t 



(A22) 



Our result can be summarized as follows. We begin 
with our total Hamiltonian, Eq. (|A1(I , which can be writ- 
ten as a positive-definite real-symmetric matrix whose 
states live in an enlarged Hilbert space containing elec- 
tromagnetic flux densities and internal oscillator modes. 
When we "integrate out" the non-resonant internal os- 
cillator modes of the media, we are left with a set of 
effective constitutive relations of the form 



ViX 



3 



B ij -(w A )u+e <UAt + c.c, 



(A23) 



and an effective Hamiltonian (which represents the con- 
served time-averaged energy density of the electromag- 
netic fields as well as the oscillator modes) that involves 
a different tensor Bij (uj) given in lA22l Using the relation 
in Ea. HA13|) . we can equivalently write the Hamiltonian 
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as 

ff=^E-By(«i)V (A24) 

For the case of generalized frequency dependence con- 
sidered here, the normalization of the electromagnetic 
fields are given (up to a scale factor) in terms of the 
time-averaged energy density, Ea. fMty : 

V (u^)*, B-^wJuJ) = —V- (A25) 

Jit ' W f 

Finally, the matrix B _1 and not B~ x enters the expres- 
sion for the Berry connection, since it also defines the 
normalization of our states. 

APPENDIX B: NUMERICAL ALGORITHMS 
FOR BANDSTRUCTURE CALCULATIONS 

In this section, we shall describe our formulation of the 
photonic bandstructure problem which has been used in 
the computations of the edge mode spectra. 

Since we always neglect absorption/emmision and 
other non-linear processes of light (i.e. we work within an 
approximation of photon number conservation), we seek 
a real-space Hamiltonian formulation of the bandstruc- 
ture problem. A real-space method is desirable over ex- 
isting Fourier-space methods for our purposes; the modes 
we are particularly interested in are obtained in domain 
wall configurations of the Faraday "mass term" as a func- 
tion of position, and it is most simple and suitable to 
work within a real space formulation. 

In the numerical implementation of a Hamiltonian for- 
mulation, we shall treat the continuum flux densities 
(y\ = (D,B) rather than (r\ = (E,H) as our funda- 
mental dynamical variables. The former set obey the 
source-free Gauss' relations: 

V • \y) = 0. (Bl) 

The Hamiltonian of our system is given by the follow- 
ing quadratic form: 

H= l -(D,e- 1 D) + 1 -{B,^B). (B2) 

Furthermore, the propagating solutions of Maxwell's 
equations require the fields to be coupled in non- 
canonical Poisson bracket relations: 

{D a (x), B b (x')} = e abc V c S 3 (x - x') (B3) 

The two sets of fields are related by \y) = B\r), where 
B is the matrix of constitutive relations introduced in 
section II. The source free Maxwell equations are slight 
variants of the ones described in section II. Written as a 
generalized eigenmode problem of form 

AB\y)=tu\y). (B4) 



The matrix A is the imaginary anti-symmetric matrix 
introduced in section II, and B = B 1 is a positive- 
definite Hcrmitian matrix. The eigenmode problem here 
is formally analogous to the problem of a non-canonical 
harmonic oscillator with Hamiltonian 

n=\Y,B ijyiyj , (B5) 
and Poisson brackets 

{y i ,y j } = A ij , (B6) 

Since A is imaginary and anti-symmetric, its eigenval- 
ues are either zero, or come in pairs with opposite sign. 
It is the presence of zero modes which prevents a canon- 
ical treatment of the problem. In the Maxwell problem, 
one third of the A matrix eigenvalues are zero modes. 




FIG. 18: The generic discretization scheme for the photon 
Band structure problem. Space is broken up into polyhedra. 
Local electric energy density is defined at the vertices of each 
polyhedron, and the electric fluxes, defined on the edges of the 
polyhedron, connect two electric energy sites. The volume of 
the polyhedron is associated with local magnetic energy den- 
sity, and magnetic fluxes "live" on the faces of the polyhedron. 
The scheme here has electric-magnetic duality in that a dual 
polyhedron can be defined on the vertices of which magnetic 
energy density defined, etc. The scheme here is inspired by 
lattice QED, which ensures the correct long wavelength pho- 
ton dispersion; the only difference here is absence of sources. 

In the spatial discretization of this problem, we di- 
vide space into polyhedral cells, whose vertices contain 
the local electrical energy density as well as the inverse 
permittivity tensor e~ 1 (r). The electrical fluxes, $£> are 
defined on the edges of the polyhedron, while the mag- 
netic fluxes, $b are associated with the faces of each 
cell. Finally, magnetic energy and the local inverse per- 
meability tensor fj,^ (r) are defined on the centers of each 
polyhedron. 

This discretization scheme preserves the self-duality of 
the source-free Maxwell equations in three dimensions; 
for each such electric polyhedron described above, there 
is a dual magnetic polyhedron whose faces correspond to 
the edges of the electric polyhedron, and whose center 
corresponds to the vertices of the electric polyhedron. 



The discretized form the A matrix couples electric 
fluxes to magnetic ones, and vice-versa. The coupling 
is (see Fig. H3| 

A? B = = 0,±i. (B7) 




FIG. 19: The discretized form of the A, which contains the 
Poisson bracket relations of the fluxes. Shown here are exam- 
ple configurations of {<&f , $ ^ } = +i (top), —i (middle), and 
(bottom). 

The B matrix couples fluxes of the same type, and de- 
pends on the geometry of the polyhedra used to discretize 
space. For the case of a simple cubic discretization, and 
for the electric fluxes (see Fig. I2U[) . 

Bu = l(eu 1 (r 1 )+e^ 1 (r 2 )) (B8) 

Ba = \ejhr2). (B9) 

Identical relations involving the inverse permeability 
tensor are constructed for the magnetic fluxes. With 
the present formulation, the complete Hamiltonian of 
the system is expressed as a sum of local terms, H — 
J2n h ( r n), with 

h(r n ) = J2 B a(rn)ViVi' (BIO) 

ij 

Using this method, we can handle the case where the con- 
stitutive relations have generalized anisotropy, and vary 
in space. 
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FIG. 20: The discretized form of the B, which contains the 
contains the geometric as well as the dynamics information. It 
couples fluxes of the same type only, and allows for anisotropy 
in the constitutive relations. 

We have made use of a simple cubic discretization of 
this general algorithm (although similar implementations 
using BCC and FCC lattices have given the correct count 
of the long-wavelength photon modes) to compute pho- 
tonic band structures. The fluxes Q D and Q B are made 
to obey the generalized Bloch boundary conditions 

<S> a (x + R) = e ikR <f>"(x),a = D,B, (Bll) 

where R is a lattice translation of the particular photonic 
crystal under consideration, and k is a Bloch vector in 
the first Brillouin zone. We compute the bandstructure of 
the system by varying the Bloch vector in the boundary 
conditions, which introduces Bloch phase factors into a 
few off-diagonal elements of the A and B matrices and 
give rise to the band dispersions. 

Both the A and B matrix are sufficiently sparse and 
are stored as matrix-vector multipliers and are treated 
using a Lanczos algorithm. However, due to the sig- 
nificantly large number of zero modes, a conventional 
Lanczos treatment of Ea. (|B4() would not converge. The 
Lanczos adaptation for the Photonic problem is done by 
modifying the A matrix to 

A -> A' = ABA - 2w A, (B12) 



A'B\y) =w(w-2uo)\y). (B13) 

Here, loq is the lowest eigenvalue, and the low lying 
(negative) eigenvalues of the modified problem are now 
the physically relevant ones which are easily found with 
the Lanczos implementation. The dimensions of the ma- 
trices are d = 6N, where N is the number of points used 
to discretize the constitutive relations. We have found 
system sizes up to 10 6 to be accessible within this ap- 
proach. Furthermore, local polarization and magnetiza- 
tion modes can be added to the algorithm. 
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APPENDIX C: DERIVATION OF THE DIRAC 
POINT SPLITTING 



The new eigenvalue problem with the symmetry break- 
ing terms is 



In this section, we derive a general expression for the 
frequency splitting at the Dirac point caused by inver- 
sion or time-reversal symmetry breaking perturbations. 
We will use "Bra-ket" notation to represent our eigen- 
vectors instead of writing equations for each component. 
We suppose that we know the exact eigenstates of the 
problem 



A\u Q ) = oj d B q 1 \u ), 



(CI) 



and that the solutions are two fold degenerate at the 
Dirac point, as for example, in the numerical examples 
we have considered. Now add a perturbation in the con- 
stitutive relations: 



B 1 =B " 1 +\B^\ 



(C2) 



This term represents our inversion or time-reversal 
breaking perturbation. To find the splitting of the Dirac 
point (our "Dirac mass"), we solve the modified problem 



A\u) = uj (Bq 1 + AJBf x ) \u). 



(C3) 



Since the B 1 matrix is positive-definite, it has a well 

defined positive-definite inverse square root matrix B\J 2 
and we can rewrite the unperturbed problem in the form 
of a conventional Hermitian eigenvalue problem 



Bl /2 ABl /2 \z ) 



v D \z ), 



(C4) 



where 



A\u) = lu (Bq 1 + A-Bj" 1 ) \u) 

= uB^ 2 (l + XB^B^B^)B^\u), 

which subsequently is rewritten in the canonical form as 



Bl /2 (A - XojB^ 1 ) Bq/ 2 \z) = uj\z), 



(C6) 



where \z) = B, 



-1/2 



u). The correction to the spectrum 



to first order in perturbation theory in the eigenvalue 
problem above is then 



-uj D \{z \Bl /2 B^Bl /2 \z ) 
(u \B^\u ) 



(U \B Q \u ) 



We have restored the normalization factor for the state 
| Mo) in the last line above. Thus, our main result here is 
a general expression for the splitting of the Dirac point 
frequency spectrum, given by the dimensionless quantity 



5uj_ = x (u \B 1 1 \u ) 
lo d (u q \Bq 1 \u ) 



(C7) 



\z ) = B- 1/2 \u ) 



(C5) 



